The Journal of Chemical Physics

Recent advances in hydrodynamic theory have revealed the severe effect of periodic boundary conditions (PBCs) on the diffusive dynamics of lipid membranes in molecular dynamics simulations. Even when accounting for PBC effects, the corrected lipid diffusion coefficients often severely overshoot the experimental estimates. Here, we investigate the underlying reasons for the exaggerated dynamics, and suggest potential ways for improvement. To this end, we examine the diffusion of four lipid types in both bilayers and monolayers using the CHARMM36 force field. We account for PBC effects using the full hydrodynamic treatment: for bilayers we use non-equilibrium simulations to extract the interleaflet friction parameter used in the correction; whereas monolayer hydrodynamics are treated by setting this parameter to zero. Our results suggest that the dynamics of bilayers are too fast, even if interleaflet friction is accounted for. However, the change of the water model to OPC leads to an excellent agreement with experiments. For monolayers, the dynamics with the TIP3P water model agree well with experiments, whereas they are undershot with OPC. As OPC and TIP3P differ in both shear viscosity and surface tension, we develop two new mass-scaled water models to clarify the roles of the thermodynamic and kinetic properties of the water model on lipid dynamics. Our results indicate that both of these quantities play a major role in lipid dynamics. Moreover, it seems that the accurate description of diffusion in both lipid bilayers and monolayers cannot be accounted for by changes in the water model alone, but likely also requires modifications in the lipid model.

Many biologically relevant processes occur on time scales that are beyond the reach of atomistic simulations. These processes include large protein dynamics and the self-assembly of biological materials. Coarse-grained molecular modeling allows computer simulations on length and time scales 2-3 orders of magnitude larger than atomistic simulations, bridging the gap between the atomistic and mesoscopic scales. However, the structural information involving the atomic planes is lost in coarse-grained. We develop a simple coarse grained protein model with structural information in explicit solvent. We represent each residues center of mass as a polymer bead and water oxygen as a solvent bead. Each polymer bead has five degrees of freedom: position of the center and additional two variables, for the backbone dihedral angles. All interaction parameters for bonded, non-bonded and dihedral coupling are derived from the equilibrated all-atom molecular dynamics simulation trajectory. We find that our coarse-grained approach, reproduces residue-level structural features that closely match the crystal structures and all-atom simulation results for both structured and disordered proteins.

Molecular dynamics simulations of mixtures of the ceramide N-(tetracosanoyl)-sphingosine (NS), cholesterol, and a free fatty acid are performed to gain a molecular-level understanding of the structure of the lipids found in the stratum corneum layer of skin. A new coarse-grained model for cholesterol, developed using the multistate iterative Boltzmann inversion method, is compatible with previously developed coarse-grained forcefields for ceramide NS, free fatty acid, and water, and validated against atomistic simulations of these lipids using the CHARMM force field. Self-assembly simulations of multilayer structures using these coarse-grained force fields are performed, revealing that a large fraction of the ceramides adopt extended conformations, which cannot occur in the bilayer structures typically studied using simulation. Cholesterol fluidizes the membrane by promoting packing defects and it is observed that an increase in cholesterol content reduces the bilayer height, due to an increase in interdigitation of the C24 lipid tails, consistent with experimental observations. Through the use of a simple reverse-mapping procedure, a self-assembled coarse-grained multilayer system is used to construct an equivalent structure with atomistic resolution. Simulations of this atomistic structure are found to closely agree with experimentally derived neutron scattering length density profiles. Significant interlayer hydrogen bonding is observed in the inner layers of the atomistic multilayer structure that are not found in the outer layers in contact with water or in equivalent bilayer structures. These results identify several significant differences in the structure and hydrogen bonding of multilayer structures as compared to the more commonly studied bilayer systems, and, as such, highlight the importance of simulating multilayer structures for more accurate comparisons with experiment. These results also provide validation of the efficacy of the coarse-grained forcefields and the framework for multiscale simulation.

The dynamics of the aggregated light-harvesting complex (LHCII) associated with its antennae pigments can be crucial for a transition between light harvesting and dissipative states pivotal for non-photochemical quenching (NPQ). To this end, aggregation of chlorophyll-a (CLA) without the LHCII and pigment binding LHCII monomers in the plant thylakoid membranes have been investigated using coarse-grained molecular dynamics simulations at 293 K. Both CLA without the LHCII and pigment-binding LHCII monomers dynamically form and break dimers and higher-order aggregates in thylakoids within the simulation time. The contact lifetime and waiting time distributions of CLA dimers exhibit multiple time scales including most populated fast time scales and less populated slow time scales. The survival probability of CLA dimer in the absence of the LHCII follows a non-exponential decay with multiple residence time scales, leading to a time-dependent rate, unlike conventional rate theory. Such non-exponential decay of survival manifests the emergence of dynamic disorder in CLA without the LHCII resulting from the coupling between time scales of dimer formation and higher-order aggregates. The conformational fluctuations of the LHCII known for inter-CLA coupling variation occur on multiple time scales comparable to the LHCII dimer residence time scales leading to less probable but comparable and more probable slower inter-CLA fluctuations. This indicates the dynamic coupling in the LHCII conformations and their aggregates with the antennae pigments can result in dynamic disorder which will be highly relevant for the light-harvesting efficiency and regulation of NPQ. TOC Graphic O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=111 SRC="FIGDIR/small/638782v1_ufig1.gif" ALT="Figure 1"> View larger version (57K): org.highwire.dtl.DTLVardef@5fd895org.highwire.dtl.DTLVardef@8456cforg.highwire.dtl.DTLVardef@5f5c69org.highwire.dtl.DTLVardef@abf6d0_HPS_FORMAT_FIGEXP M_FIG C_FIG

Biomolecular condensates are important contributors to the internal organization of the cell material. While initially described as liquid-like droplets, the term biomolecular condensates is now used to describe a diversity of condensed phase assemblies with material properties extending from low to high viscous liquids, gels, and even glasses. Because the material properties of condensates are determined by the intrinsic behaviour of their molecules, characterising such properties is integral to rationalising the molecular mechanisms that dictate their functions and roles in health and disease. Here, we apply and compare three distinct computational methods to measure the viscoelasticity of biomolecular condensates in molecular simulations. These methods are the shear stress relaxation modulus integration (SSRMI), the oscillatory shear (OS) technique, and the bead tracking (BT) method. We find that, although all of these methods provide consistent results for the viscosity of the condensates, the SSRMI and OS techniques outperform the BT method in terms of computational efficiency and statistical uncertainty. We, thus, apply the SSRMI and OS techniques for a set of 12 different protein/RNA systems using a sequence-dependent high-resolution coarse-grained model. Our results reveal a strong correlation between condensate viscosity and density, as well as with protein/RNA length and the number of stickers vs. spacers in the amino-acid protein sequence. Moreover, we couple the SSRMI and the OS technique to nonequilibrium molecular dynamics simulations that mimic the progressive liquid-to-gel transition of protein condensates due to the accumulation of inter-protein {beta}-sheets. We compare the behaviour of three different protein condensates--i.e., those formed by either hnRNPA1, FUS, or TDP-43 proteins--whose liquid-to-gel transitions are associated with the onset of amyotrophic lateral sclerosis and frontotemporal dementia. We find that both SSRMI and OS techniques successfully predict the transition from functional liquid-like behaviour to kinetically arrested states once the network of inter-protein {beta}-sheets has percolated through the condensates. Overall, our work provides a comparison of different modelling rheological techniques to assess the viscosity of biomolecular condensates, a critical magnitude that provides information on the behaviour of biomolecules inside condensates.

Biological events occurring on long timescales, such as protein folding, remain hard to capture with conventional molecular dynamics (MD) simulation. To overcome these limitations, enhanced sampling techniques can be used to sample regions of the free energy landscape separated by high energy barriers, thereby allowing to observe these rare events. However, many of these techniques require a priori knowledge of the appropriate reaction coordinates (RCs) that describe the process of interest. In recent years, Artificial Intelligence (AI) models have emerged as promising approaches to accelerate rare event sampling. However, integration of these AI methods with MD for automated learning of improved RCs is not trivial, particularly when working with undersampled trajectories and highly complex systems. In this study, we employed the State Predictive Information Bottleneck (SPIB) neural network, coupled with bias exchange metadynamics simulations (BE-metaD), to investigate the unfolding process of two proteins, chignolin and villin. By utilizing the high-dimensional RCs learned from SPIB even with poor training data, BE-metaD simulations dramatically accelerate the sampling of the unfolding process for both proteins. In addition, we compare different RCs and find that the careful selection of RCs is crucial to substantially speed up the sampling of rare events. Thus, this approach, leveraging the power of AI and enhanced sampling techniques, holds great promise for advancing our understanding of complex biological processes occurring on long timescales. O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=139 SRC="FIGDIR/small/550401v1_ufig1.gif" ALT="Figure 1"> View larger version (29K): org.highwire.dtl.DTLVardef@18a9296org.highwire.dtl.DTLVardef@9dc2f3org.highwire.dtl.DTLVardef@16a0cf8org.highwire.dtl.DTLVardef@1799db1_HPS_FORMAT_FIGEXP M_FIG TABLE OF CONTENT GRAPHIC C_FIG

Brownian dynamics (BD) simulations that include hydrodynamic interactions (HIs) modeled at the Rotne-Prager-Yamakawa (RPY) level of theory are a valuable tool for accurately modeling the translational and rotational diffusion of macromolecules such as proteins and nucleic acids. A major drawback to the inclusion of HIs in BD simulations is their computational expense, and an obvious way to consider reducing the expense of BD-HI simulations is to include a cutoff such that HIs beyond a certain distance are omitted. Unfortunately, a naive attempt to implement such a scheme usually leads to the RPY diffusion tensor becoming non-positive definite, which has the consequence that it becomes impossible to compute the correlated random displacements required by the Ermak-McCammon BD-HI algorithm. Here I show that a simple approach can be used to overcome this problem and implement a distance-based cutoff scheme that is guaranteed to lead to a diffusion tensor that is positive definite. The method involves only a straightforward distance-based scaling of the original RPY terms, and allows a seamless transition to be made between BD simulations that neglect HIs entirely and simulations that include HIs at the full RPY level of theory.

We provide a molecular-level description of the thermodynamics and mechanistic aspects of drug permeation through the cell membrane. As a case study, we considered the anti-malaria, FDA approved drug chloroquine. Molecular dynamics simulations of the molecule (in its neutral and protonated form) were performed in the presence of different lipid bilayers, with the aim of uncovering key aspects of the permeation process, a fundamental step for drugs action. Free energy values obtained by well-tempered metadynamics simulations suggest that the neutral form is the only permeating protomer, consistent with experimental data. H-bond interactions of the drug with water molecules and membrane headgroups play a crucial role for permeation. The presence of the transmembrane potential, investigated here for the first time in a drug permeation study, does not qualitatively affect these conclusions. TOC Graphic O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=200 SRC="FIGDIR/small/550356v1_ufig1.gif" ALT="Figure 1"> View larger version (71K): org.highwire.dtl.DTLVardef@d1f5b9org.highwire.dtl.DTLVardef@5b6b0forg.highwire.dtl.DTLVardef@1d1b55borg.highwire.dtl.DTLVardef@d96b52_HPS_FORMAT_FIGEXP M_FIG C_FIG

1Stochasticity due to fluctuations in chemical reactions can play important roles in cellular network-driven processes. Although the Stochastic Simulation Algorithm (SSA, aka Gillespie Algorithm) has long been accepted as a suitable method to solve the time-dependent chemical master equation, its computational cost is prohibitive for large scale complex networks such as those found in cellular processes. Here we present GPU-SSA, an implementation of the SSA formalism utilizing Graphics Processing Units for use in Python using the PySB modeling framework. We show that the GPU implementation of SSA can achieve significant speedup compared to parallel CPU or single-core CPU implementations. We further include supplementary didactic material to demonstrate how to incorporate GPU-SSA workflows for interested readers.

Covalently attached sugar molecules play important roles as mediators of biomolecular interactions. Molecular dynamics simulations are an indispensable tool to explore these interactions at the molecular level. The large time and length scales involved frequently necessitate the use of coarse-grained representations, which heavily depend on the parameterization of sugar-protein interactions. Here, we adjust the sugar-protein interactions in the widely used Martini 2.2 force field to reproduce the experimental second virial coefficients between sugars and proteins. In simulations of two model proteins in glucose solutions with adjusted force field parameters, we observe weak protein-sugar interaction. The sugar molecules are thus acting mainly as crowding agents, in agreement with experimental measurements. The procedure to fine-tune sugar-protein interactions is generally applicable and could prove useful also for atomistic force fields.

The pathways by which a molecular process transitions to a target state are highly sought-after as direct views of a transition mechanism. While great strides have been made in the physics-based simulation of such pathways, the analysis of these pathways can be a major challenge due to their diversity and variable lengths. Here we present the LPATH Python tool, which implements a semi-automated method for linguistics-assisted clustering of pathways into distinct classes (or routes). This method involves three steps: 1) discretizing the configurational space into key states, 2) extracting a text-string sequence of key visited states for each pathway, and 3) pairwise matching of pathways based on a text-string similarity score. To circumvent the prohibitive memory requirements of the first step, we have implemented a general two-stage method for clustering conformational states that exploits machine learning. LPATH is primarily designed for use with the WESTPA software for weighted ensemble simulations; however, the tool can also be applied to conventional simulations. As demonstrated for the C7eq to C7ax conformational transition of alanine dipeptide, LPATH provides physically reasonable classes of pathways and corresponding probabilities. TOC Graphic O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=111 SRC="FIGDIR/small/553774v2_ufig1.gif" ALT="Figure 1"> View larger version (15K): org.highwire.dtl.DTLVardef@14eed4corg.highwire.dtl.DTLVardef@bd1f67org.highwire.dtl.DTLVardef@58c04borg.highwire.dtl.DTLVardef@b88034_HPS_FORMAT_FIGEXP M_FIG C_FIG

Computational modeling and simulation of biomolecular systems at their functional pH ranges requires an accurate approach to exploring the pH dependence of conformations and interactions. Here we present a new approach - the Equilibrium Constant pH (ECpH) method - to perform conformational sampling of protein systems in the framework of molecular dynamics simulations in an N, P, T-thermodynamic ensemble. The performance of ECpH is illustrated for two proteins with experimentally determined conformational responses to pH change: the small globular water-soluble bovine b-lactoglobulin (BBL), and the dimer transmembrane antiporter CLC-ec1 Cl-/H+. We show that with computational speeds comparable to equivalent canonical MD simulations we performed, the ECpH trajectories reproduce accurately the pH-dependent conformational changes observed experimentally in these two protein systems, some of which were not seen in the corresponding canonical MD simulations. O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=153 SRC="FIGDIR/small/394015v1_ufig1.gif" ALT="Figure 1"> View larger version (79K): org.highwire.dtl.DTLVardef@b9c6bdorg.highwire.dtl.DTLVardef@a6d0d3org.highwire.dtl.DTLVardef@1d7fb23org.highwire.dtl.DTLVardef@a8f760_HPS_FORMAT_FIGEXP M_FIG Table of Contents artwork C_FIG

Using molecular dynamic simulations we study whether amyloidogenic regions in viral proteins can initiate and modulate formation of -synuclein aggregates, thought to be the disease-causing agent in Parkinsons Disease. As an example we choose the nine-residue fragment SFYVYSRVK (SK9), located on the C-terminal of the Envelope protein of SARS-COV-2. We probe how the presence of SK9 affects the conformational ensemble of -synuclein monomers and the stability of two resolved fibril polymorphs. We find that the viral protein fragment SK9 may alter -synuclein amyloid formation by shifting the ensemble toward aggregation-prone and preferentially rod-like fibril seeding conformations. However, SK9 has only little effect of the stability of pre-existing or newly-formed fibrils.

Since it was first observed, the COVID-19 pandemic has created a global emergency for national health systems due to millions of confirmed cases and hundreds of thousands of deaths. At a molecular level, the bottleneck for the infection is the binding of the receptor binding domain (RBD) of the viral spike protein to ACE2, an enzyme exposed on human cell membranes. Several experimental structures of the ACE2:RBD complex have been made available, however they offer only a static description of the arrangements of the molecules in either the free or bound states. In order to gain a dynamic description of the binding process that is key to infection, we use molecular simulations with a coarse grained model of the RBD and ACE2. We find that binding occurs in an all-or-none way, without intermediates, and that even in the bound state, the RBD exhibits a considerably dynamic behaviour. From short equilibrium simulations started in the unbound state we provide snapshots that result in a tentative mechanism of binding. Our findings may be important for the development of drug discovery strategies that target the RBD.

Weak magnetic field exposure can affect many biological processes across a wide range of living organisms. Recently, it has been observed that weak magnetic fields can modulate reactive oxygen species (ROS) concentration, affecting regeneration in planaria. These effects show unusual nonlinear dependence on magnetic field strength, including a sign change. In another study by the same group, superoxide is identified as the particular ROS being modulated. We propose a radical pair mechanism based on a flavin-superoxide radical pair to explain the modulation of superoxide production and its effect on planarian regeneration. The results of our calculations favor a triplet-born radical pair. Our yield calculations can reproduce the observed magnetic field dependence, including the sign change. Moreover, to explain the size of the effect on ROS concentration, we suggest a simple amplification model inspired by known biochemical mechanisms and lay out the conditions for such a model to work. Further, we also make empirical predictions concerning the hypomagnetic field effects on planarian regeneration.

An efficient variational method is presented for estimating the diffusion coefficients and free-energy profiles along selected collective variables from projected molecular dynamics trajectories under both equilibrium and nonequilibrium conditions. The method is based on the assumption that the short-time transition probability density of the coordinate moves can be approximated by a Gaussian form. Defining a loss function as the sum of Kullback-Leibler divergences between the analytical short-time propagators of an overdamped Langevin model and those estimated directly from the projected trajectories maximises the agreement between the two and allows for its analytic evaluation. To efficiently minimise this loss function by varying diffusion and free-energy profiles along collective variables, we use an adaptive Monte Carlo scheme. The method is applied to two model systems exhibiting diffusive dynamics, as well as to water diffusion across the interface of a biomolecular condensate, demonstrating its robustness and accuracy.

We perform all-atom molecular dynamics simulations to study the effects of the N-linked glycans on the stability of the spike glycoprotein in SARS-CoV-2. After a 100 ns of simulation on the spike proteins without and with the N-linked glycans, we found that the presence of glycans increases the local stability in their vicinity; even though their effect on the full structure is negligible. Graphical Abstract O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=95 SRC="FIGDIR/small/475397v1_ufig1.gif" ALT="Figure 1"> View larger version (32K): org.highwire.dtl.DTLVardef@cb1c1dorg.highwire.dtl.DTLVardef@a2ecbeorg.highwire.dtl.DTLVardef@64db68org.highwire.dtl.DTLVardef@180a1aa_HPS_FORMAT_FIGEXP M_FIG C_FIG

The effects of Guanidine Hydrochloride (GdmCl) on two Intrinsically Disordered Proteins (IDPs) are investigated using simulations of the Self-Organized Polymer-IDP (SOP-IDP) model. The impact of GdmCl is taken into account using the Molecular Transfer Model(MTM). We show that, due to dramatic reduction in the stiffness of the highly charged Prothymosin- (ProT) with increasing concentration of GdmCl ([GdmCl]), the radius of gyration (Rg) decreases sharply till about 1.0M. Above 1.0M, ProT expands, caused by the swelling effect of GdmCl. In contrast, Rg of -Synuclein (Syn) swells as continuously as [GdmCl] increases, with most of the expansion occurring at concentrations less than 0.2M. Strikingly, the amplitude of the Small Angle X-ray Scattering (SAXS) profiles for ProT increases till [GdmCl]{approx} 1.0M and decreases beyond 1.0M. The [GdmCl]-dependent SAXS profiles for Syn, which has a pronounced bump at small wave vector (q [~] 0.5nm-1) at low [GdmCl] ([≤] 0.2M), monotonically decrease at all values of [GdmCl]. The contrasting behavior predicted by the combination of MTM and SOP-IDP simulations may be qualitatively understood by modeling ProT as a strongly charged polyelectrolyte with nearly uniform density of charges along the chain contour and Syn as a nearly neutral polymer, except near the C-terminus where the uncompensated negatively charged residues are located. The precise predictions for the SAXS profiles as a function of [GdmCl] can be readily tested.

Lipid domains in cellular membranes can adopt either registered or anti-registered configurations between the two leaflets, playing crucial roles in numerous cellular processes. However, the microscopic mechanisms governing domain registration and anti-registration remain incompletely understood due to the complexity of interleaflet interactions. In this work, we investigate the effects of lipid headgroup size and cholesterol concentration on domain registration using coarse-grained molecular dynamics simulations. Through a free energy perturbation (FEP)-like method, we systematically demonstrate that headgroup size and cholesterol concentration cooperatively regulate domain registration by modulating the membrane curvature and line tension at domain boundaries. Furthermore, cholesterol flip-flop strengthens interleaflet coupling, thereby facilitating domain registration. This work quantitatively reveals how lipid geometric properties and cholesterol concentration and translocation jointly regulate interleaflet coupling, offering novel insights into the molecular mechanisms underlying membrane domain organization. TOC Graphic O_FIG O_LINKSMALLFIG WIDTH=200 HEIGHT=114 SRC="FIGDIR/small/699701v1_ufig1.gif" ALT="Figure 1"> View larger version (43K): org.highwire.dtl.DTLVardef@ef52b7org.highwire.dtl.DTLVardef@1fa026org.highwire.dtl.DTLVardef@4b2d36org.highwire.dtl.DTLVardef@185659d_HPS_FORMAT_FIGEXP M_FIG C_FIG

The need to incorporate specific molecular-scale features for largescale structural changes in biological membranes necessitate use of a multi-scale computational approach. Here, this comprises of Langevin dynamics in a normal mode space determined from an elastic network model (ENM) representation for lipid-water Hamiltonian. All atom (AA) MD simulations are used to determine model parameters, and Langevin dynamics predictions for an extensive set of bilayer properties, such as, undulation spectra, undulation relaxation rates, dynamic structure factor, and mechanical properties are validated against the data from MD simulations and experiments. The transferability of model parameters to describe dynamics of a larger lipid bilayer and a heterogeneous membrane-protein system is assessed. The developed model is coupled to the energy landscape for membrane deformations to obtain a set of generic reaction coordinates (RCs) for pore formation in two tensionless, single lipid-type bilayers, namely, 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) and 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC). Structure evolution is carried in an AA MD simulation wherein the generic RCs are used in a path metadynamics or an umbrella sampling simulation to investigate thermodynamics of pore formation and its molecular determinants. The transition state is characterized extensively to bring out the interplay between various bilayer motions (undulations, lateral density fluctuations, thinning, lipid tilt), lipid solvation, and lipid packing.